Predicting High-Density Polyethylene Melt Rheology Using a Multimode Tube Model Derived Using Non-Equilibrium Thermodynamics

Based on the Generalized bracket, or Beris–Edwards, formalism of non-equilibrium thermodynamics, we recently proposed a new differential constitutive model for the rheological study of entangled polymer melts and solutions. It amended the shortcomings of a previous model that was too strict regarding the values of the convective constraint release parameter for the model not to violate the second law of thermodynamics, and it has been shown capable of predicting a transient stress undershoot (following the overshoot) at high shear rates. In this study, we wish to further examine this model’s capability to predict the rheological response of industrial polymer systems by extending it to its multiple-mode version. The comparison with industrial rheological data (High-Density Polyethylene resins), which was based on comparison with experimental data available in (a) Small Amplitude Oscillatory shear, (b) start-up shear, and (c) start-up uniaxial elongation, was noted to be good.


Introduction
As reported by the Society of Plastics Industries (SPI) in 2000, the plastic industry in the US is positioned, in terms of shipment, in fourth place among manufacturing industries, following motor vehicles and equipment, electronic components and accessories, and petroleum refining [1]. A more recent survey predicts that the global plastic packaging market will be worth $269.6 billion by 2025, achieving a compound annual growth rate (CAGR) [2]. This figure alone highlights the impact of plastic materials on our lives and, thus, the significance of optimizing polymer processing technology. Future polymer processing will focus not on the machine, but on the product [1]. Several instabilities appear in the polymer industry, which makes life difficult for polymer engineers. For example, under certain circumstances, when molten plastic is forced through a die, the shark-skin defect appears [3]. To avoid this issue, it is suggested to slow down the manufacturing rate; however, this action decreases the production rates of commercial products, leading to an increase in cost. Wang et al. have suggested that this defect may be related to a molecular instability that corresponds to an oscillation of the absorbed chains in the die exit area between coiled and stretched states [4]. Thus, it seems that the answer needed should be sought by maintaining the molecular level of description and performing molecular dynamics simulations. The ultimate goal of this process is to predict the properties of a product via numerical simulations based on molecular principles and multiple-scale techniques [1].
Due to computational limitations, however, this goal has been unachievable until the last few years, in which period the extended evolution of simulation algorithms, the parallelization of these algorithms and the race, which has been very recently undertaken, that both models performed equally well (note that the thermodynamic admissibility of these two models is shown in Ref. [32]). Hoyle et al. [33] evaluated the performance of the multimode pom-pom model to those of both LDPE and branched HDPE melts, whereas, more recently, Konaganti et al. [15] employed the double convected pom-pom [34] model to predict the rheological behavior of a high-molecular-weight HDPE melt. Since multimode versions have been found to be superior to single-mode versions, our aim in this paper is to generalize the tube model of Stephanou et al. [23] to its multimode version and use it to predict the rheological response of a HDPE melt.
The structure of the paper is as follows: In Section 2, the new model is briefly derived using non-equilibrium thermodynamics (NET). In Section 3, we derive the expressions of the relevant rheological material functions obtained by analyzing the asymptotic behavior of the model in the limits of small shear rates. The results obtained with the new model are then presented in Section 4, in which we first discuss its parameterization and then show how accurately and reliably it can describe the viscoelasticity of HDPE polymer melts. The paper concludes with Section 5, in which the most important aspects of our work are summarized, and future plans are highlighted and discussed.

State Variables
This work considers an isothermal and incompressible flow, meaning that the total mass density ρ and the entropy density (or temperature) are excluded from the vector of state variables. To characterize the polymer chains, the entanglement strand conformation tensor c, following the method of Stephanou et al. [23,35], is considered to be made dimensionless through c = Kc/k B T, where K denotes the spring constant of the Hookean dumbbells that represents the entanglement strands at equilibrium, k B the Boltzmann constant, and T denotes the absolute temperature [36]. The conformation tensor c refers to one entanglement strand, and at equilibrium (zero flow field applied), it coincides with the unit tensor. To characterize the multiple modes of the polymer chains, N conformation tensors are considered, with one tensor being considered for each mode [36]. Finally, the momentum density M, which is the hydrodynamic variable, is further considered, meaning that, overall, the vector of state variables is expressed as x = {M, c (1) , c (2) , . . . , c (i) , . . . , c (N) }.

System Hamiltonian
The mechanical part of the system's Hamiltonian is given as [36] where represents the kinetic energy of the system, whereas [23,35,36] represents the system's Helmholtz free energy (with a(x) the Helmholtz free energy density) that includes the following contributions: the dimensionless potential Φ tr c (i) − I , which accounts for chain stretching, and an entropic contribution, which involves the logarithm of the determinant of the conformation tensor of each mode [35,36]. Here, G (i) e is the entanglement modulus of the i th mode, and I is the unit tensor. The partial derivative of the potential with respect to the trace of the conformation tensor defines the (dimensionless) effective spring constant [37,38] for the i th mode meaning that the corresponding Volterra derivative of the free energy with respect to the conformation tensor is Here, the following FENE-P(Wagner) expression is used: e is the finite extensibility (FENE) parameter of the entanglement strand associated with the i th mode. As shown by the study of Stephanou et al. [35], b e = 3 (0.82) 2 /C ∞ (M e /M 0 ) (when all FENE parameters are considered equal), where C ∞ is the polymer characteristic ratio at infinite chain length, M e is the entanglement molecular weight, and M 0 is the average molar mass of a monomer. For example, for PS melts, b e = 54 [35].

The Poisson and Dissipation Brackets
Following the work of Beris and Edwards [36], the Poisson bracket for multiple modes is given as follows: We note both here and throughout this paper that Einstein's summation convention for repeated Greek indices is employed. The complete Poisson bracket was then simply given as follows: Next, the following expression for the dissipation bracket associated with the conformation tensors is used [36].
The first integral on the right-hand side of Equation (6) accounts for the relaxation effects of each conformation tensor, which is proportional to a fourth-rank relaxation tensor, whereas the second integral allows for the non-affine deformation of each conformation tensor. We note that the subscript "nec", meaning "no entropy production correction", is added to the dissipation bracket to indicate that this bracket lacks terms that involve Volterra derivatives with respect to the entropy density, which can be omitted when we consider isothermal systems [36]. We further note that although, in general, the dissipation bracket allows explicit coupling between cross modes, as shown in Equations (8.2-25) of Beris and Edwards [36], in this study, they are omitted for simplicity. Then, where we have defined the upper-convected time derivative: .
Finally, the extra (polymeric) stress tensor is given as follows:

The Matrices L and Λ
The relaxation matrix of each mode Λ (ii) αβγε is split into two contributions, following the work of Stephanou et al. [23], which have different relaxation times: Here, τ (i) CR is the CR relaxation time of the i th mode, which is considered to be half of the corresponding reptation time, τ (i) [39] [we note that this time coincides with the CCR relaxation time at equilibrium, as shown in Equation (10)], and τ (i) R (tr c) is the Rouse relaxation time of the i th mode, R,eq is the equilibrium Rouse relaxation time of the i th mode, which is given as τ [19], and k (i) is the Extended White-Metzner (EWM) exponent [36] for the i th mode. We note that for the Rouse time, a shear rate dependency through the use of the trace of the conformation tensor of each mode is considered. The functions f Rouse tr c (i) are scalar functions of the trace of the conformation tensor, as defined via the following equation [23]: ccr is the CCR parameter of the i th mode. For the (dimensionless) mobility tensor β (i) of the i th mode, the Giesekus' postulate β (i) = I + α (i) σ (i) is used [37] with σ (i) = σ (i) /G i , and α (i) is the anisotropic mobility (Giesekus) parameter of the i th mode. Then, with Λ where [23,39] the CCR relaxation time is obtained as follows: Finally, the expression of the L (i) matrix is given via the following equation [36,37]: Here, ξ (i) is the non-affine/slip parameter of the i th mode. This parameter is important, as it allows for the prediction of a transient stress undershoot (following the overshoot) at high shear rates [23].

Thermodynamic Admissibility
Any thermodynamic system must obey the restriction of a non-negative rate of total entropy production. When the fluid studied is isothermal and incompressible, the entropy production results from the degradation of the mechanical energy, leading to [36]. For this aspect to be satisfied in our model, the following equation must hold: ., N} are the three eigenvalues of the conformation tensor of the i th mode. Obviously, since the conditions 0 ≤ α [23] guarantee that each term of the summation is positive, the sum as a whole is also positive, meaning that the multimode version of the Stephanou et al. model [23] presented in this work is thermodynamically admissible.

Conformation Tensor Evolution Equation
The evolution equation used for each of the dimensionless conformation tensors is obtained by substituting Equations (4b), (9d), and (11) The CCR relaxation time for the i th mode is given in Equation (10), and the (dimensionless) effective spring constant is given in Equation (4c). Finally, the expression for the polymeric stress tensor is obtained by substituting Equations (4b) and (11) into Equation (8) as follows:

Asymptotic Behavior of the Model in Steady State Shear
In this section, we provide analytical expressions that describe the asymptotic behavior of the multimode version of the Stephanou et al. [23] model in the limit of weak flows for the following cases: inception of simple shear flow (SSF) described by the kinematics u = γ cos(ωt)y, 0, 0 , in which ω is the oscillation frequency. The material functions to be analyzed are as follows: (a) the transient shear viscosity η ε), in the case of uniaxial elongation; and (c) the storage, G (ω), and loss, G (ω), moduli in the case of SAOS.
To obtain the asymptotic behavior, we need to expand the conformation tensor for each mode in the limit of small strain rates (by invoking a linearization of the evolution equation for each of the conformation tensors) and analytically solve the corresponding ordinary differential equations. After this stage, we obtain the non-zero stress tensor components via Equation (14). Finally, we obtain the following results for the relevant material functions: Inception of shear: Inception of uniaxial elongation: meaning that Trouton's law is true for the steady-state extensional viscosity. Small Amplitude Oscillatory Shear: In Equations (15)-(17), we have defined G

Results and Discussion
The FENE parameter can be easily calculated via b e = 3(0.82) 2 M e /(M 0 C ∞ ), as mentioned above. In this study, PE M 0 = 14 g/mol, whereas C ∞ = 7.3 and M e = 828 g/mol (see Table 3.3, p. 151 of Ref. [40]). These values yield b e = 16.34. We will compare against the experimental data of Konaganti et al. [15] that have performed rheological measurements of the sample HDPE-1 (reported by the same group [41]), for which M w = 206 kg/mol; thus, the number of entanglements is equal to Z ≈ 249 >> 1. The relaxation spectrum is the same as the one used by Konaganti et al. [15] (see their paper's Table 2 for T = 200 • C, though is also provided in Table 1), and it was obtained by fitting the expressions of the storage and loss moduli, which are shown in Equation (17), with the corresponding experimental data. The comparison against the experimental storage and loss moduli is shown in Figure 1. . Figure 1. Comparison between the model predictions Eqs. (17) and the experimental data presented in Ref. [15] for the storage and loss moduli of an HDPE sample at 200 °C. The relaxation spectrum is provided in Table 1.
All of the remaining parameter values are obtained by fitting the model predictions with the experimental data; for simplicity, we assume that each parameter has the same value for all modes, although, in general, different values for each mode could be considered (e.g., Konaganti et al. [15]). The following values of the model parameter are chosen to provide a good comparison with the start-up shear flow experimental data provided in Figure 2: ξ = 0.02, α = 0.3, βccr = 4 × 10 −4 , and k = −3.5. For comparison, we also depict, in the following figures, the predictions of the multimode Giesekus model [which is a special case of our model in which βccr = ξ = k = 0 and be infinite or the function h = 1 in Equation (4c)] with α = 0.3 and the relaxation spectrum of Table 1. Figure 2(a) illustrates the comparison between the experimental data for the growth of the shear viscosity upon inception of shear flow and the simulated results obtained using the new model. The experimental data (symbols) were collected at three different shear rates-0.05 1/s, 0.5 1/s, and 1 1/s-whereas the lines represent the simulated shear viscosity values at the corresponding shear rates. It can be observed that the model accurately captures the trends and magnitude of the shear viscosity over time, doing so much more successfully than the Giesekus model [ Figure 2(b)]. We should, however, note that the overshoots noted at the two larger shear rates (0.5 1/s and 1 1/s) are higher than the experimental data. The overshoot predicted using our model is controlled by two parameters: the FENE parameter be and the anisotropic mobility (Giesekus) parameter α. As mentioned above, the former parameter is not a free parameter, as it is directly dictated by structural parameters. The latter parameter is a free parameter, and by increasing its value, the start-up shear viscosity overshoots noted at the two larger shear rates (0.5 1/s and 1 1/s) shift downwards and broaden (results not shown), thus more closely agreeing with the experimental data; however, the good comparison identified at the smaller shear rate (0.05 1/s) is reduced. This result might hint that the parameter α should not be a constant, but should increase with the applied strain rate. Similar arguments have been put forth and resulted in a variable non-affine/slip parameter proposed by Nikiforidis et al. [42] and a variable link tension coefficient proposed by Stephanou and Kröger [26]. We note that although a non-zero value of ξ is employed, the undershoots produced are too small to be noted via the scale used in Figure 2. Although no experimental data are provided, we provide the corresponding prediction of the growth of the first and second Figure 1. Comparison between the model predictions Equation (17) and the experimental data presented in Ref. [15] for the storage and loss moduli of an HDPE sample at 200 • C. The relaxation spectrum is provided in Table 1. All of the remaining parameter values are obtained by fitting the model predictions with the experimental data; for simplicity, we assume that each parameter has the same value for all modes, although, in general, different values for each mode could be considered (e.g., Konaganti et al. [15]). The following values of the model parameter are chosen to provide a good comparison with the start-up shear flow experimental data provided in Figure 2: ξ = 0.02, α = 0.3, β ccr = 4 × 10 −4 , and k = −3.5. For comparison, we also depict, in the following figures, the predictions of the multimode Giesekus model [which is a special case of our model in which β ccr = ξ = k = 0 and b e infinite or the function h = 1 in Equation (4c)] with α = 0.3 and the relaxation spectrum of Table 1. Figure 2a illustrates the comparison between the experimental data for the growth of the shear viscosity upon inception of shear flow and the simulated results obtained using the new model. The experimental data (symbols) were collected at three different shear rates-0.05 1/s, 0.5 1/s, and 1 1/s-whereas the lines represent the simulated shear viscosity values at the corresponding shear rates. It can be observed that the model accurately captures the trends and magnitude of the shear viscosity over time, doing so much more successfully than the Giesekus model (Figure 2b). We should, however, note that the overshoots noted at the two larger shear rates (0.5 1/s and 1 1/s) are higher than the experimental data. The overshoot predicted using our model is controlled by two parameters: the FENE parameter b e and the anisotropic mobility (Giesekus) parameter α.

Comparison with Start-Up Shear Flow Data
As mentioned above, the former parameter is not a free parameter, as it is directly dictated by structural parameters. The latter parameter is a free parameter, and by increasing its value, the start-up shear viscosity overshoots noted at the two larger shear rates (0.5 1/s and 1 1/s) shift downwards and broaden (results not shown), thus more closely agreeing with the experimental data; however, the good comparison identified at the smaller shear rate (0.05 1/s) is reduced. This result might hint that the parameter α should not be a constant, but should increase with the applied strain rate. Similar arguments have been put forth and resulted in a variable non-affine/slip parameter proposed by Nikiforidis et al. [42] and a variable link tension coefficient proposed by Stephanou and Kröger [26]. We note that although a non-zero value of ξ is employed, the undershoots produced are too small to be noted via the scale used in Figure 2. Although no experimental data are provided, we provide the corresponding prediction of the growth of the first and second normal stress coefficients in Figure 3, as well as the steady-state values of all viscometric material functions in Figure 4.         Figure 2. The multimode Giesekus model prediction is also depicted (black dashed lines).

Comparison with Start-Up Uniaxial Elongational Flow Data
In Figure 5, we provide the comparison between the model predictions and the experimental measurements obtained by Konaganti et al. [15] for the growth of the elongational viscosity as a function of time upon inception of uniaxial elongation at the following stretch rates: 0.05 1/s, 0.5 1/s, and 5 1/s. The comparison employs the same parameter values as those used in Figure 2, except ξ = 0, as informed by Stephanou et al. [37], and the corresponding steady-state prediction is provided in Figure 6. Given that the parameter values were selected based on the start-up shear data (Figure 2), the comparison is noted to be adequately appropriate. We note that the Giesekus model (Figure 5b) fails to reproduce correctly the trend of the experimental data. It should be noted that the uniaxial elongation data obtained by Konaganti et al. [15] do not seem to reach a steady state. This outcome is customary in the literature [11,12,43] due to experimental difficulties, as the polymer samples used become too slender and often break. This issue has led some researchers to consider the highest value as the steady-state elongational viscosity and omit all following data [44]. Thus, the data that follow the overshoot may not be accurate. Overall, the proposed model demonstrates good agreement with the experimental data, indicating its effectiveness in describing the rheological behavior of polymers of industrial significance.

Comparison with Start-Up Uniaxial Elongational Flow Data
In Figure 5, we provide the comparison between the model predictions and the experimental measurements obtained by Konaganti et al. [15] for the growth of the elongational viscosity as a function of time upon inception of uniaxial elongation at the following stretch rates: 0.05 1/s, 0.5 1/s, and 5 1/s. The comparison employs the same parameter values as those used in Figure 2, except ξ = 0, as informed by Stephanou et al. [37], and the corresponding steady-state prediction is provided in Figure 6. Given that the parameter values were selected based on the start-up shear data (Figure 2), the comparison is noted to be adequately appropriate. We note that the Giesekus model ( Figure 5b) fails to reproduce correctly the trend of the experimental data. It should be noted that the uniaxial elongation data obtained by Konaganti et al. [15] do not seem to reach a steady state. This outcome is customary in the literature [11,12,43] due to experimental difficulties, as the polymer samples used become too slender and often break. This issue has led some researchers to consider the highest value as the steady-state elongational viscosity and omit all following data [44]. Thus, the data that follow the overshoot may not be accurate. Overall, the proposed model demonstrates good agreement with the experimental data, indicating its effectiveness in describing the rheological behavior of polymers of industrial significance.  [15] for the growth of the elongational viscosity as a function of time for several stretch rates. The thick black line depicts the LVE envelope, as given in Equation (16).
The sparameter values are the same as those used in Figure 2 except ξ = 0. In panel (b), we depict the same comparison for the multimode Giesekus model.

Conclusions
In this work, we developed the multiple-mode version of a generalized, conformation-tensor based viscoelastic model for polymer melts, which was proposed in [23], by making use of the generalized bracket formalism of Beris and Edwards [36]. Like its forerunner, i.e., the single-mode version [23], it accounts for the most significant effects that can be realized in entangled polymer systems: anisotropic drag, finite extensibility, non-affine motion (leading to the exhibition of a transient stress undershoot at large shear rates), and variable chain relaxation due to convective constraint release. The multiplemode version of the model has been shown to have a very good predictive capability with regard to the industrial experimental data for HDPE obtained by Konaganti et al. [15]. Obviously, a better comparison could have been obtained if we had simultaneously fitted Figure 5. (a) Comparison between the model predictions (lines) and the experimental measurements (symbols) of Konaganti et al. [15] for the growth of the elongational viscosity as a function of time for several stretch rates. The thick black line depicts the LVE envelope, as given in Equation (16). The sparameter values are the same as those used in Figure 2 except ξ = 0. In panel (b), we depict the same comparison for the multimode Giesekus model.  [15] for the growth of the elongational viscosity as a function of time for several stretch rates. The thick black line depicts the LVE envelope, as given in Equation (16).
The sparameter values are the same as those used in Figure 2 except ξ = 0. In panel (b), we depict the same comparison for the multimode Giesekus model.

Conclusions
In this work, we developed the multiple-mode version of a generalized, conformation-tensor based viscoelastic model for polymer melts, which was proposed in [23], by making use of the generalized bracket formalism of Beris and Edwards [36]. Like its forerunner, i.e., the single-mode version [23], it accounts for the most significant effects that can be realized in entangled polymer systems: anisotropic drag, finite extensibility, non-affine motion (leading to the exhibition of a transient stress undershoot at large shear rates), and variable chain relaxation due to convective constraint release. The multiplemode version of the model has been shown to have a very good predictive capability with regard to the industrial experimental data for HDPE obtained by Konaganti et al. [15]. Obviously, a better comparison could have been obtained if we had simultaneously fitted

Conclusions
In this work, we developed the multiple-mode version of a generalized, conformationtensor based viscoelastic model for polymer melts, which was proposed in [23], by making use of the generalized bracket formalism of Beris and Edwards [36]. Like its forerunner, i.e., the single-mode version [23], it accounts for the most significant effects that can be realized in entangled polymer systems: anisotropic drag, finite extensibility, non-affine motion (leading to the exhibition of a transient stress undershoot at large shear rates), and variable chain relaxation due to convective constraint release. The multiple-mode version of the model has been shown to have a very good predictive capability with regard to the industrial experimental data for HDPE obtained by Konaganti et al. [15]. Obviously, a better comparison could have been obtained if we had simultaneously fitted all available data, as Konaganti et al. [15] did (we note that some parameters must have different values in different flows, such as the non-affine/slip parameter, which must be explicitly considered in the fitting process). Furthermore, different values of the model parameters for each mode could also have been considered, following the work of Konaganti et al. [15], which would certainly provide more flexibility to the model and, thus, improve its capacity to fit the experimental data. However, in our present work, we mostly focused on deriving the multimode version of the model proposed by Stephanou et al. [23] using non-equilibrium thermodynamics, with less focus devoted to its predictive capacity to almost perfectly fit all available experimental data.
The model, in its present form, only considers strictly linear chains. Industrial samples, particularly LDPE, are never strictly linear, having either short or long branches distributed along their entire backbone. There is clear evidence that all material functions of PE are considerably affected by the presence of even low levels of long-chain branching [45]. We, therefore, need to generalize it and allow the explicit consideration of branches, following the guidelines provided by the pom-pom [27] model and its thermodynamically admissible version, known as the pom-pon [28] model. Furthermore, the multimode version does not explicitly consider the molecular weight (MW) distribution of industrial samples, such as the log-normal or gamma distributions, that are able to describe experimentally measured distributions [46]. As such, we should also generalize it to handle molecular weight distribution, following the work of Schieber [47,48]. This generalized constitutive model will allow a more accurate prediction of the rheological responses of industrially used polymeric systems that possess an extensive spectrum of MW. Finally, we used an ambiguous model-fitting process wherein the values of the model parameters were obtained to best compare them to the experimental data. However, atomistic non-equilibrium molecular dynamics simulations can be used to directly obtain some of the parameters from the simulations [37,49]. For example, the EWM exponent can be obtained by directly comparing the prediction of Equation (9b) to the relaxation time as a function of shear rate at large shear rates obtained via the NEMD simulations, and the CCR parameter can be obtained by directly comparing the average number of entanglements as a function of shear rate which, due to CCR, is noted to decrease [22] (we note, however, that in our present work we omitted this approach and considered a constant number of entanglements). Only then, polymer engineers would be able to accurately use the predictions of the multimode constitutive model for comparison with the rheological response noted in actual industrial processes. Our findings provide a foundation for future research that aims to enhance the properties of high MW polymers for diverse applications.